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EQUATIONS FOR DETERMINING AIRCRAFT MOTIONS FROM ACCIDENT DATA 
Ralph E. Bach, Jr., and Rodney C. Wingrove 
Ameo Research Center 


SUMMARY 


This paper describes procedures for determining a comprehensive accident 
scenario from a limited data set. The analysis techniques accept and process 
data from either of two sources: an Air Traffic Control (ATC) radar-tracking 

system or a foil flight-data recorder. Local meteorological information at 
the time of the accident and aircraft performance data are also utilized. 

Equations for the desired aircraft motions (roll, pitch, yaw) and forces 
(lift, excess thrust) are given in terms of elements of the measurement set 
and certain of their time derivatives. The principal assumption made in the 
development presented here is that aircraft side force and sideslip angle 
are negligible. 

An estimation procedure is outlined for use with each data source. For 
the foil case, a discussion of exploiting measurement redundancy is given. 
Since either formulation requires estimates of measurement time derivatives, 
an algorithm for least-squares smoothing is provided 


INTRODUCTION 


For several years Ames Research Center has been assisting the National 
Transportation Safety Board (NTSB) and the military services in their investi- 
gations of aircraft accidents. During this period, methods have been developed 
to determine aircraft motions from the limited data available following an 
accident. It is the purpose of this report to derive the equations used in 
the analysis of such data and to summarize their application. A companion 
paper (ref. 1) presents an experimental assessment of the accuracy of the 
methods and results from several accident analyses. To the authors' knowledge, 
these are the first expositions of accident analysis techniques based on the 
full aircraft kinematic model that have appeared in the open literature. 


The data sources considered here are the ground-based radar tracking 
system and the onboard foil flight recorder. A radar tracking system provides 
time histories of aircraft position, including altitude (when the vehicle is 
equipped with an altitude transponder). Radar data are recorded by an Air 
Traffic Control (ATC) station in the vicinity of the accident site (ref. 2). 

A foil record contains time histories of indicated airspeed, magnetic heading, 
barometric altitude, and normal acceleration. Flight recorders of the metal- 
foil type are carried by many airliners and some military aircraft. They are 
of rugged construction, designed to withstand the rigors of severe accidents 
(refs 3 and 4). 
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By 1974 , NTSB scientists had developed digital computer programs for 
analysis of both ATC radar data and foil records# these programs corrected 
the data for calibration errors and nonstandard atmospheric conditions and 
provided time~hlstory estimates of such parameters as grounds peed and ground^ 
track and flightpath angles (ref# 5)# The NASA contribution has been to show 
how either data source, in combination with meteorological and aircraft data, 
can be utilised to estimate wind-axis specific forces (excess thrust, lift) 
and body-axis Euler angles (roll, pitch, yaw). Thrust and lift forces are 
Important aircraft performance parameters (ref, 6), which, along with the 
Euler angle time histories, are valuable aids to the investigator in visual- 
ising an accident trajectory. 

This report is organised as follows} 1) wind-axis force and angle rela<-. 
tions are derived for a vehicle governed by basic rigid-body dynamics, 
referred to a flat, nonrotating earth} 1) then, formulas for expressing body- 
axis Euler angles in terms- of wind-axis angles and tl e aerodynamic angles are- 
presented; 3) the specialised applications of the equations to the analysis 
of either radar or foil data are then summarised; pertinent axis systems and 
transformations are defined in Appendix A; and finally, since both analyses 
depend on estimates of first and second derivatives of certain measurement 
time histories, an algorithm for least-squares smoothing is presented in 
Appendix B.. 


FORCE ANDUATTITUDE-EQUATIONS 


The key steps in processing accident data from ATC stations or foil 
recorders are indicated in figure 1. Smoothing the data provides estimates 
of first and second time derivatives fro m w hich wind— axis forces and angles 
Can be determined, The forces, airspeed, and specific information about the 
airplane (lift characteristics) permit estimating the angle of attack, which 
can then be used with the wind-axis angles to determine time histories of the 
body-axis Euler angles. The assumption of negligible side force (ayw - 0) 
and sideslip angle (0 ■ 0) is necessary to obtain the desired solutions. 

There may be, of course, accident situations in which this assumption is not 
valid. In the following paragraphs, we derive the equations for force and 
attitude estimation. The vehicle dynamics Introduced here are discussed in 
reference 7. 


Wind-Axis Forces and Angles 


It is convenient first to develop expressions for wind-axis angles and 
forces in terms of quantities derived from vehicle position derivatives and 
local wind estimates. The rationale for this approach should become evident 
as the development proceeds. Orientation of the wind-axis frame F w , relative 
to the vehicle-carried vertical frame Fy, is defined by Euler angles (<|» w , 6 W , 
$ w ). As described in Appendix A, both Fy and Fy have their origins at the 
vehicle mass center and move with it. Frame Fy has axis Oxy along the 
velocity vector and axis Ozy in the vehicle plane of symmetry. The axes of 
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Figure 1.— Flow chart for processing accident data to obtain wind- 
axis forces and body-axis Euler angles. 
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Fy are parallel to the stationary earth-surface frame Fg. In this develop- 
ment we assume that the earth's curvature and rotation relative to an Inertial 
frame can be neglected. Definitions of the Euler- angles and a review of 
Important axis transformations can also be found In Appendix A. 

Expressions for wind-axis angles and 0 W are easily derived from the 
relationship between vehicle velocity in Fg and Fy, given by 


( 1 ) 


where (k t £, i) ere time derivatives of vehicle position- in E., (w v . w«. w*) 
are winds in F E , and V is true airspeed. It should be noted that, along 
a space trajectory, winds would be ideally characterized.. as functions of time 
and position, i.e., w* - Wx<t,.x,y,z); w y - w y (t,x,y,z); w z - w z (t,x,y,z). 

When the second equation of (1) is divided by the £irst, it is seen that 


1 
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V cos0 w Cosi|> w 
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■ hnt 
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V cos6w slntJv 
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0 
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-V sin6 w 


tartt|» w - (y - w y )/(x - w x ) 

(2) 

while the third equation yields 


sin6 w - -(£ - w z )/V 

(3) 


Note that the airspeed can be written as 


v " [(* “ w x ) 2 + - w y ) 2 + (& - w z ) 2 ] 1/2 (, 

which is obtained by premultiplying the vector in equation (1) by its trans- 
pose. 


The relationship for vehicle aerodynamic reaction in terms of the ac- 
celeration of its center of mass and its gravitational attraction is used to 
complete the wind-axis fiuler angle set and determine the forces* Hence, 


a xw 
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Aerodynamic reaction as defined in reference 7 is 


a xw 


T xw “ D 

a yw 

- 

^yw “ c 

a zw 

m 


T ZW - L 

Ml * 


(5) 


( 6 ) 


where (Txw, Ty*,, T zw ) are thrust terms; D is drag; C is side force; L is 
lift} and m is vehicle mass. For an aircraft in cruising flight, only 
T xw> D > and L are significant forces. Often, is referred to as "excess 

thrust." The acceleration of the vehicle center of mass can be written 
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a cw ■ r wV *y O'. 

-q w V 

where $ and (tf x , ifry, ti B ) are time derivatives of airopeed and winds* The 
angular velocities q w and r w are given by 

q w " 4 w ain K + K coa0 w cofl ^w ( 
r w " § w C08 *w + *w UOfl0 w 8in *w 1 

Note that expressions for the derivatives 4» w , $ w , V can be obtained by 
differentiating equations (2) through (4), which yields 

^w " " w y )costji w - (x - * x )sint|> w J/V cose w (9] 

6 w « - (z - <* z + V sinO w )/V cos0 w (10; 

V - [«- w x )eosiJ> w + (y - & y )sini|i w ]co80 w - (2 - * z )sinO w (11] 

where (x, y, z) are second time derivatives of vehicle position in Fg. 

We now multiply both sides of equation (5) by L x (-'|> w ) and use equations 
(7) » (8), and (10) to obtain 

■*, • ♦ + <8 - W 2 )sin0 w + (w x cosi|> w + w y sim|» w )cose w (12) 

and 

-*w 8ln *w + V COB *w ■ C 1 | 

V 8ln *» + *»» C08 *w ■ C 2 

where 

C 1 “ V coar w + (w y cos * w - w x 8imj» w ) ^ 

C 2 - (2 - g + a^ sin0 w )/cos0 w | 

Next, we solve uquation (13) for sln<J> w , cos<J> w as 

8l "*w ■ «V„ - C X“ Z W )/( V + “£> { (U) 

co 8 * w - (C l0jrw * C 2 a w >/(a£ + n z 2 w ) | 

from which we obtain the exproBBion 


m ( a yw^ a xw) ^ (C1/-C2) 

'w * 1 - (a /a Xcj-Cn) 
yw zw 7 x 1 
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It should be noted that side reaction flyw usually small compared with 
a ?w . In fact, In order to estimate 4> w and a 0W by using the formulation of 
equations (13) and (14), it is necessary to assume that the ratio (ny^n^ 
is negligible. In that case 


tan<|> w ■ (C^/-C 2 ) 

Finally, aolutlon of equation (13) for a gw yields 

«« ' °2 ““♦w - C 1 9ln *w 

An analysis of the orrors in iw and a aw incurred by noglocting the ratio 
(Qy^/ a z ^) shows that 


( 17 ) 


(18) 


tan^ • 

(a c /a M ) - 1 - ,«♦, 


! 


(19) 


where <|>e l 8 error bbe estimate of <|> w and 
estimate of a zw , 


a e is the error in the 


Body-Axis Euler Angles 

We have shown that wind-axis Euler angles (1> w > 9w» +w) and specific 
forces (a^, a w , a zw ) may be determined from vehicle position derivatives 
Ind wind 2 forSiioE It is, however the attitude of the body frame Fg 
relative to F v as described by the Euler angles (<|q>, 0b* ♦b) tha * ia of 
primary interest in the analysis of vehicle motion. Recall that the frames 
L and y F W are displaced by angle of attack o and sideslip angle 8 (see 
Appendix A). In this development, we assume that sideslip is negligible, 
which is consistent with the statement a yw * 0. In the next section, we 

obtain an estimate of a by using a xw> a zw> V and X should be 

particular aircraft. Hence, given (Vw* °w* *w) aad 5° To 

possible to determine expressions for the body-axis angles O^b* b> ^b'* 
conclude this section, therefore, we shall review the transformation betwe 
frames Fg and Fy and derive the expressions relating body and wind-axis 

Euler angles. 


The body-axis system Fg is carried into the wind-axis system Fw by 
the rotation sequence (0, -a, 0). The transformation matrix is given by 



Hra “ 

L z (0)L y (-oO 

(20) 

Now, note that a vector 

b in Fg 

may be transformed into 

Fy in either of 

two equivalent ways: 

v " Wb 15 

or v » LyyLyj, b 

(2L) 


6 
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\ 
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^jsisSaitss^BM 


Hare v in n vector with coordinate* in Fy »nd Wb repreaente the trnnn- 
formation 

Nra " V-VV-W"*^ <22> 


Hone© wo can ertprona 


Wb " HwHfB 


(23) 


n s , md / a g\ wo (Should bo able to solve for (4>b* 0b* 
J”) in 8 racr to'do U ;i.l", wo ten She matrix product indlontod by o,untlon 

(22) t 


t m 


^Vb 


COS0. cos«(». 

ain<l> b sinG b coaiJ' b 

co8<(> b 8in0 b coB<J» b “ 

D D 

- cos<j> b sin^ b 

+ 8in<(> b sin^ b 

cos0. siwiv 

sin$ b 8in9 b situ|> b 
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D D 

•f C08(|> b COB<|» b 

- sin<)) b cosi|> b 

-sin0 u 

» D 

sin<t> b coa0 b 

COS<j> b CO80 b 


(24) 


from which it follows that 


tan<t> b 

sin0 b 

tani|< b 


*32^33 

“ a 31 

*2i/ A u 


(25) 

(26) 
(27) 


‘X - E> «“» *‘'" ati ° n <27) to obtatn cha “ pr<i 


tan (4*^ - V - <*a n°d*„ - Hi 8ltl *w )/U ll c ° e *w + *21 8l °*w> 


(28) 


Finally, we evaluate the mat rix prod uct of ^“ at J°“ 0 (2 ^ a ^one°™5)” (M)! ‘ 
£ (SrXtS: formulae 


tan<|> 


cos0 din<j' w cosO^ - sin 8 sln®^ 


b - - elnn ainB elnd s )C 0 Bb M - einn coaB eln» M 

alnb b -coac cobB aln8 w + (etna coBd w + cobc einB Bln* M )coB8 w 

sino sin<t> w “ cosa sinp cobi}> w 

tan(t|» b -1^)- c os« co7e ~^r~ - (sina cob4> w + cosa sing sin^sinT 


(29) 

(30) 
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(31) 


TABLE X, - ELEMENTS OF MATRIX I,, 
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MOTIONS DERIVED PROM ACCIDENT DATA 


In this section, procedures^ for determining aircraft motions and forces 
from radar or foil data arc sunmafiaeds the procedures are baaed on the 
equationo developed in the proviouo section, We nnaumo that certain typea of 
data anomalies, such aa dropouts, wild points, etc., havo been removed, and 
that corrections for temperature, proaoure, etc,, havo been made to air data 
in no far an possible. Meteorological data, of course, will t • impreci.-e - 
wind measurements, if available, may havo boon recorded many milea from the 
crash site and at a different time from that of tho accident. Usually, data 
sufficient to determine only wx and Wy (an functions of altitude) are 
obtained. Hence, it must often bo assumed that w 2 “ 0, These limitations 
notwithstanding, we can then proceed to compute wind-axis Euler angles and 
forces along the trajectory, to estimate an angle-of-uttack-time history, and 
to determine the body-axis Euler angles. The necessary assumption for a 
solution with either data sot is that vehicle side force fly W and sideslip 
angle E are negligible, 


Analysis of Radar Data 

We assume that time hi? Tories of vehicle position (x, y, z) nr * mds 
(wx* Wy, w 2 ) are available. First, "smooth" the data (see 4i>tn noi . h) to 
obtain estimates of vehicle velocity (x, y, i ) and acceleration (K. y, E). 

We then use equations (2) through (4) to determine wind-axis angles iji w and 
0 W and airspeed Vs 

>J> W “ tan -1 [(y - w y )/(x - w x >] 

6 w “ sin~^£-(li - w z )/v] 

v - [(*. - w x ) 2 + - w y r + <* - w z ) 2 ] 1;/2 

Next, we determine excess thrust by using equation (11) in equation (12), 
which yields 


\ 


(32) 


i w - (« cosi|> w + V slnt|» w )cose w - (E - g)sin0 w 

anA then we solve for lift and wind-axis roll angle with 

V ■ C 2 C08 *w - C 1 sin *v 1 
♦ w " ton“ 1 (C 1 /-C^) J 

where equation (9) has been used in equation (14) to give 

- 9 coss|> w - a 8ir.* w i 

C 2 - (E - g + a^ Sin0 w )/cos9 w ) 


(33) 


(34) 


( 35 ) 


Notice that in equations (33) and (35) there is no dependence on wind- 
derivative (fc x , & y , $%) information. 

Finally, In order to obtain the body-axis angles (i)^, fljj, 4>b)» on estimate 
of angle of attack is needed. In the linear region, that estimate can bo 
approximated by 


a » a 0 - m a^/QS Cl 
Q - (1/2 )pV 2 


(36) 


where m Is mass, S is wing area, p Is air density, and CL a Is the 
derivative of the lift coefficient with respect to a. Values of a 6 and C» 
depend primarily on aircraft configuration and Mach number and are tabulated 0 
for a given aircraft. For high angles of attack, a flat-plate relationship 
yields a good approximation (refs.— $L-and 10): 

a - tan“ 1 (a xw /a 2w ) (37) 


Now the body angles can be determined by using equations (29) through (31), 
with 6 ■ 0: 

<b ■ + tan -1 jsina (cOsot cosO^ - sina cos^ sln0 w )J \ 

0 b « sin""‘*‘(cosa sin0 w + sina cos0 w cos<f> w ) > (38) 

<b. » tan”'*' [cos0 sin<|> / (cosa cos0„ cos<|> - sina sin0 )"] / 

b Lww ww wj 


Analysis of Foil Data 

The data set derived from a foil record consists of true airspeed V, 
heading angle altitude -z, and normal specific force ~ a %h’ These 

quantities, along with time derivatives V, <J»b> z, and z, wind information, 
and aircraft performance data, can be utilized to provide an accident scenario. 
First, we use equation (3) to determine the wind-axis pitch angle 0 W : 

© w • sin ^jj-(fe - w z )/vj (39) 

Next, the forces a^ and a 2w and the roll angle <fr w are estimated by using 
equations (12) through (14) : • 

fljjw “ ^ + (g - w z )sin0 w + (& x cos* w + sin^ w )cos0 w (40) 

e *v " c 2 COO K ~ c i Bin K 

<P W - tan“ 1 (C 1 /-C 2 ) 


( 41 ) 


C t - if' w V W1H0 W + * y ™«4' w - * x Hln^ | (42) 

U 2 m (« - H + <» xw HlnO w )/aoHO w 

Notice that the exprinwIonH for a xw , n zw , ami >|> w depend ou wliuKU>rlvat I ve 
quantities (* x , * v , V* «'»tlroatoH of which are honor a. I ly unreliable (or 
unavailable) . 8 taco w,, in usually not measured, ami w x , w y are reeouled 

aa func t Iona only of altitude, tho boat estimates ol the. wind dorlvat Ivoa an 

w a * ft; \sf a ft 5 W„ " II ^ , l) 

w x dir. * y «« « 

Not Leo also that the factor <J' W in equation (42) muat bo approximated by the 
time derivative of f> h , whloh La included In tho data act:. ( I wind-dot’ l vat Ive 
offoota In equations (40) and (42) are to bo Included, lot <J- W " % alao.l . 
Finally , after estimating anp.lo of attack aa In the radar caae, wo ho Ive for 

the unknown Kulor any, lea: 


0 ® aln” 1 (o.oau alnO w + slnu coaO w cos«fr w ) 

+ - tmr l [coaO w altti|» w / (coaa co«0 w coh* w - alna HlnO w >] 


w'J i 

coh4* _)1 * 


1'w " ^b “ tan" l [aln« ain«|> w / (coao eoaO w - ainu «lnO w eos* w )] * 

The foil record apparently contains redundancy, alnoo tho measurement 
a.., has not appeared in any of the oquatlona. The redundancy may be utilized 
to provide an independent estimate of angle of attack. In the linear region 
that estimate, can be approximated by 

« - % - ma 2b /QSC Na 

where m, Q, and S have been previously defined and is the derivative 

of the normal force coefficient. Again, o» n and would be found as 

tabulated function*) of configuration and Mach number for a given aircraft. 
Outside of the linear region, a may be obtained from a flat-plate relation 
equivalent to equation (37) , which is 

a » sitt~ l <V /a zb ) (Ab> 

Using the estimate of « so obtained, we can determine the wind-axis reaction 
term a zw independently of That expression is given by 

» zu - «,b /c08 “ - v tam (47) 

One way that tho redundancy Implied by equation (47) can be utlllacd 1» tn 
provide a check on the assumption that the ratio (a yw /a zw ) is negligibly sma . 
If each of the relations in equation (15) Is squared and added, the result l 


2 2 

V + “w 


C l + C 2 
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from whleh the map.nl tilde 




+ 




•> 

#- 

i'.W 



In obtained, Tho redundnnoy ran alao I**' uned to ooroponuaio tor a noma 1 1 on lu 
l ho data. Onboard lnnt rumen ta tiro nubjoet to uuuaual opera t Inp, eondltlona 
dut l up an aoo Idont , and port. Iona of tho rooovorod lol l rooord may contain 
alp, nil loant orrora. II tho alt (motor la nuapoet , t.hon tho rol l angle i|» w 
ran bo oomputod from 

♦ w * a In" l l~C,/a,, w ) < r,0 > 


When data from tho dl root tonal gyro are conaldoPod unrol lablo, 
ot d> w ran bo determined from tho oxproauLau- 




*] 


1/2 


tho magnitude 
(SI) 


Tho alternate aolutlona for roll angle pi von by equationa (SO) and (St) 
oon.ld . ot oourao, bo uaod to provide a data-oona latonoy ohook of tho four foil 
moaanromonta , In thin regard, one further obaervatlon ahould bo noted: 
Manipulation of equation (41) ylelda 

0., « a eoa>|' (S2) 

2 i’.w 1 w 

from wh I eh an exproaaion for vert leal aeeeleratlon 

« ” OOM ^w l * OM °w “ a xw Hln0 w 4 *' (V ’ ) 


follow**. lienee, equation (ST) may bo integrated twlee, and with appropriate 
oho lee of eonatanta the re.au.lt ahould mat eh the altitude rooord. 


OONCUUHNC RKMAKKS 


Thin report: preaenta the oquatlona neeoaaary to eouatruet a eomprohenalvo 
aeenario of alreraft mot Iona by ualnp, data from an Air Truffle Control radar 
ayatem or a foil flight-data reeordor, Vohlelo per forma no e data and loon l 
met eorol ogival information are alao ut l lined. 

Kxpreaa to:u» for anglea and foreea are flrat obtained In a w l ml ax I a 
frame, In terma of the moaaured varlablea and eertaln of their time dorlvatlvea. 
Kat | mat ton of ang i» of attaek then permit a determination of tho deal rod body 
axla Kulor anp.loa. The dor l vat Iona require an anaumpttou ot negligible aide 
j or ei* and aldealtp: error oat Imatoa have boon t tie hided. 

The proooiluroa deaorlbod heroin have boon Imp l omen tod to aaalat the 
Nat ional Tranaport at Ion Satotv Hoard and the military aervloea In their In 
veat Ip.at Iona id aeveral alroralt aooldeuta. 
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AWKNIUX A 


AXltf 


HYHTHMH AND TKANSKOKMAT U'NS 


In the text, die tour reference Crimen l’V,« l ‘V » >''w» l '» au> Mtl U/m1 * 

They are defined art follows! 

la) the rth-aurfaro frame t». "I* «£• 

„!u. point inn north and Ogyg point ln« on«l . 

,„> The vehicle-carried vortical frame F V; 

oi v Vt . 

(c) Tin- wiml- axis frame P M , dimit^nlnng 

l ow, lies in tile vehicle plane of symmetry. 

W> *» body-axis frame F», 

fixed at the masa center and has its ^J t B Uie) , A xtrt 

zr^z i xrz 

in flight dynamics, the Ruler "^^ 0 ^ at^^^ 

X - 7 b rJ!?;, X «Sci^. 

i““ deseribc^tho sequence of rotations Uluntrateed 

in figure 2. 

(a) A rotation is made about Ot v , carrying the axes to °*2W- 
This ia the yaw angle 4>. 

(b) A rotation la made about Oy 2 , carrying the nxaa to 0% m z V 
This is the pitch angle 0. 

(c) A rotation 1. made about Oxj. carrying the axes to their final 
K) position Oxyz. This is the roll angle +. 
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FLIGHT PATH 



V z 3 _ 


Figure 2.- - Euler angles. 


or 


The matrices required to transform a vector from F v into either F, 
*B correspond to the sequence of rotations (i|>, 0, 4 .) and are given b' 


= W L y ( VW 


for the wind axes, and 


L BV " W L y ( VW 


for the body axes. The transformations associated with a single rotation 
about each ot the coordinate axes are 



1 

0 

0 


0 

COS(f» 

sinij> 


0 

"Sitvji 

COSt|) 


eosO 

0 

-sitid 

V ’> ' 

0 

1 

0 


sitiO 

0 

0080 
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(Al) 


(A2) 


(A3) 


(A4) 




CQflty 


I^ity) * 


-sinty 

0 


flinty 0 
entity 0 
0 1 

m 


(A5) 


it should be noted that the products L x (ty)L x ('-ty) , Ly(0)L y (-0) , and L^ifO^C-ty) 
each reduce to the unit matrix. Thus, the transformations of equations (Ai) 
and (A2) are orthogonal, that is: 


Ij VW " H-^w^^w^x^w^ " HrV 

l vb " V-VV“V Lx< “V " L BV 


where the superscript T denotes the nuxtrix tr&tispose* 

Finally, it should be noted that frames Fy and Fg are displaced by 
angle of attack a and sideslip angle 3. Specifically, the body-axis 
system Fg is carried into the wind-axis system F w by the rotation sequence 
(0, -a, (3). The transformation matrix is given by 

«*» - L 2 <B)t, y (-a) <A7) 
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APPENDIX B 


A SMOOTHING ALGORITHM 


SMOOTH is a routine that "smooths" a data record and provides estimates 
of its first and second derivatives with rospdct to the independent variable. 
The algorithm, adapted from reference 11, passes a "least-squares moving arc" 
through the data. The arc is a Second-degree polynomial Spanning an Interval 
of NS (odd) equally spaced data points such that 


+ ^(i - p) + 02(1 - p) 2 
p » j - (NS - l)/2 


(Bl) 


where (Sq, a^, a^) are chosen to minimize 

i*J 

J - 2 (2 i ' V 2 

i=j-NS+l 


(B2) 


The data record consists of the samples z., j ■ 1, NPTS. Values of the 

polynomial and its first two derivatives are computed at the central point p 
of each interval as the arc is passed through the data. Estimates near the 
beginning and end of the data record are obtained by extrapolation of the 
first and last polynomials. 


A FORTRAN IV listing of the program is given in figure 3. The subroutine 
call statement is CALL SMOOTH (Z, X, Y, W, NS, H, NPTS), where the parameters 
are 


Z 

X 

Y 

W 

NS 

H 

NPTS 


input vector of length NPTS containing data record to be smoothed 

output vector of length NPTS containing smoothed data record 

output vector of length NPTS containing first derivative of X 

output vector of length NPTS containing second derivative of X 

number of points desired in smoothing interval; must be odd and 

3 < NS < NPTS 


sample interval, sec 

number of points in data record 
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SUBROUTINE SMOOTH f X, V , w, NS# H, IMPT8) 

OIMEHSIQH KnP!8IiX(nPT81 t (NPTS) , * (imPTS) 

C 

C RfcF | NHITfc SANUS HM (IF OR MfcTmHiS, MR UO 
C . 

IF (mS#I»T , NRTS ) NSFNPTS 
C3«?,/<H**g) 

c _ .. 

C COEFFICIENT CALCULATION 
C 

Nr NS 

EN8M*n**2 

01bn*(ENSQ»«,) 

D2iEnSQ»J* - 

Ctl8«7b*(3 ( *ENSQ»7 t )/l)l 
Cl 3s-i«> # /Dl 

Cg 2«U>/ tN»D2) „ 

C33sMg # *CtS/D2 

c 

C INITIAUIE .SU M. LA LC ULAli Q N 

C 

APa;ll ) 

feP.80.1 _ 

CPan, 

C 

C COMPETE S UNS RE CURSI VE! y FOK Nb&, NS 

C ~ ' 

00 1 00 Na3# NS# 2 

— nnjjlLN-lUzi _ 

NPlSlNfl )/g 
NMlSU«NMji**8 

...NPlSUfNM l *±l . _ _ 

NXTl*N-t 

NXT?*N 

-APlAMlg ( NXTU ♦ 2 ( N XT 2 ) 

BP80M-AP*NMU7CNXTn+NP1 **(MXT2> 

1 CP8CP-g*FBP*AM'»(NMiSri)*Z (NXTl)f (NPt8l?)8Z(NXT2) 

L0J>_ JCflHJjN_U£ 

“ NSTP8NPT8»NM1 

C 

C s M0OiN|j>... ¥ A l, U ts F ,CLM_Lr U2# nP 1 

C 

AORCHiAPfCl JACP 

1 _AL«Cil.tSP.._ 


Figure 3.— Listing of smoothing routine. 
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a ' * ™ T~““ ' ' "" ” :r ' ' ' ' ’ 

in iiwib i i * m i m i Kw 



Ag«ClliAP*CII«CP 
DO 190 Jfl'NPl 
L*I-NPl 

X<n«AOH*(AWi.«A|.) 

V(t)«(A4AlaU«AI)/N 
W(I)sCl»A2 
190 CONTINUE 

lF(N|,Ea,NPTI) 00 TO BfO 

c 

C INTERIOR IHOOTHINO <!,OT,NPM I.LE.NIfP) 

C 

MfNPttl 

00 100 IIMiNITP 

l8TPI«NPt 

NXTPItNNl 

AP«AP*IU»IT)AI<NXT) 
BPtBP«AP4NHUEUST)ANPl«{(NXT) 
CP«CP*t|*9P*AP«NMlOQAiaiT)4NP180*2(NXT) 
AQ*Cli*AP«ClStCP 
AI»C>a*BP 
AgiCl J*APfC19*CP 
X(n«AO 
VCX)PAt/H 
M(l)aCl*Ag 
200 CONTINUE 
290 CONTINUE 
C 

C SMOOTHED VALUED POP I«N|fP*i,***#NPf| 

C 

MaNlTPti 

DO 100 S»M#NPTI 

L»I*NETP 

X (I)*AO+lt(Al+L«Ag) 

V(I)9(AU2 t *L«A2)PH 
NCI)tC9*Ag 
100 CONTINUE 
RETURN 
END 


Figure 3.— Concluded. 
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